slope <- map_dfr(rigal_trends,
  ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)

slope_df <- slope %>%
  pivot_wider(names_from = "response", values_from = "linear_slope")

1 Summary

summary_slope <- slope %>%
  group_by(response) %>%
  summarise(summ = list(enframe(summary_distribution(linear_slope, na.rm = TRUE)))) %>%
  unnest(cols = summ) %>%
  pivot_wider(names_from = "name", values_from = "value")
summary_slope %>%
  mutate(response = get_var_replacement()[response]) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
  scroll_box(width = "100%", height = "400px")

response	min	1st_quart	median	2nd_quart	max	mean	sd	n
Appearance	-0.0559441	-0.0018271	0.0009691	0.0065140	0.0600000	0.0027393	0.0098864	2678
Chao (binary, similarity)	-0.0969697	-0.0061754	0.0000000	0.0000000	0.1210526	-0.0029650	0.0138050	2678
Chao evenness	-1.0762433	-0.0197209	0.0000000	0.0214427	0.8698342	0.0000118	0.0687473	2678
Chao species richness	-1.0804000	-0.0363727	0.0002926	0.0473422	0.9908720	0.0073944	0.1387268	2678
Chao shannon	-0.7782364	-0.0270535	0.0007389	0.0362457	0.7822804	0.0043759	0.0822483	2678
Chao simpson	-0.6747701	-0.0241523	0.0007189	0.0304959	0.5440252	0.0028965	0.0691378	2678
Disappearance	-0.0486529	-0.0005100	0.0000000	0.0042135	0.0557576	0.0017718	0.0071868	2678
Evenness	-0.0830806	-0.0060849	0.0000000	0.0069975	0.0817728	0.0006221	0.0144629	2678
SER_a (rel abundance)	-0.1059755	-0.0146337	-0.0032189	0.0007927	0.1169447	-0.0072135	0.0183244	2678
Horn (binary, similarity)	-0.0764433	-0.0091583	-0.0027316	0.0012988	0.0859619	-0.0041715	0.0113065	2678
Jaccard (binary, similarity)	-0.0781734	-0.0115266	-0.0037510	0.0014597	0.0672796	-0.0051970	0.0135836	2678
Jaccard (binary, dissimilarity)	-0.0672796	-0.0014597	0.0037510	0.0115266	0.0781734	0.0051970	0.0135836	2678
Log species richness	-20.6068500	-0.9402512	0.0498002	1.4080463	26.9911617	0.2947554	2.6535099	2678
Log total abundance	-31.2811462	-2.8068822	0.4993106	4.1780859	55.8343888	0.8272072	6.6756922	2678
Nestedness (jaccard)	-0.0830556	-0.0041642	0.0002881	0.0074839	0.1156349	0.0021108	0.0137269	2678
Shannon	-0.2007711	-0.0088851	0.0004650	0.0121768	0.1574411	0.0016449	0.0221596	2678
Simpson	-0.0829149	-0.0045179	0.0001071	0.0058510	0.0800458	0.0007642	0.0111222	2678
Species richness	-1.6727273	-0.0328546	0.0020269	0.0536209	1.5727273	0.0148247	0.1422215	2678
Total turnover (codyn)	-0.0559441	-0.0013253	0.0032257	0.0097964	0.0621302	0.0045111	0.0115567	2678
Total abundance	-755.8336538	-1.8770744	0.1129827	2.2456287	398.4499936	1.7786740	31.5332673	2678
Turnover (jaccard)	-0.1223776	0.0000000	0.0000000	0.0077107	0.0989899	0.0030863	0.0134089	2678

bs <- map(var_temporal_trends, ~summary_distribution(slope_df[[.x]]))
names(bs) <- var_temporal_trends

Jaccard trends: -0.005 (0.014 s.d.) (-0.01 in @dornelas_assemblage_2014).
- meaning -0.05 % change per decade in average
- Jaccard decreased in 65% of the sites (79% in @dornelas_assemblage_2014)
Richness trends:
- +0.3% richness change per year (2.7 s.d.), i.e. 3% per decade.

2 Check trends

p <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
    ggplot(aes(x = linear_slope)) +
    geom_histogram() +
    labs(title = .y)
)
plot_grid(plotlist = p)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

no NA in slopes since fix rigal_classification:

na_siteid_slope <- unique(slope[is.na(slope$linear_slope), ]$siteid)

sd_num <- analysis_dataset %>%
  group_by(siteid) %>%
  summarise(across(where(is.numeric), ~(sd(.x))))

sd_num %>%
  filter(siteid %in% na_siteid_slope) %>%
  summary(.)
#>     siteid            species_nb  total_abundance log_total_abundance
#>  Length:0           Min.   : NA   Min.   : NA     Min.   : NA        
#>  Class :character   1st Qu.: NA   1st Qu.: NA     1st Qu.: NA        
#>  Mode  :character   Median : NA   Median : NA     Median : NA        
#>                     Mean   :NaN   Mean   :NaN     Mean   :NaN        
#>                     3rd Qu.: NA   3rd Qu.: NA     3rd Qu.: NA        
#>                     Max.   : NA   Max.   : NA     Max.   : NA        
#>  log_species_nb      year        latitude     longitude      main_bas  
#>  Min.   : NA    Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA  
#>  1st Qu.: NA    1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Median : NA    Median : NA   Median : NA   Median : NA   Median : NA  
#>  Mean   :NaN    Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA    3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA    Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA  
#>  chao_richness  chao_shannon  chao_simpson chao_evenness   hillebrand 
#>  Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA  
#>  1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Median : NA   Median : NA   Median : NA   Median : NA   Median : NA  
#>  Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA  
#>      total       appearance  disappearance    jaccard         horn    
#>  Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA  
#>  1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Median : NA   Median : NA   Median : NA   Median : NA   Median : NA  
#>  Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA  
#>       chao     total_abundance_int    shannon       simpson     inv_simpson 
#>  Min.   : NA   Min.   : NA         Min.   : NA   Min.   : NA   Min.   : NA  
#>  1st Qu.: NA   1st Qu.: NA         1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Median : NA   Median : NA         Median : NA   Median : NA   Median : NA  
#>  Mean   :NaN   Mean   :NaN         Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA   3rd Qu.: NA         3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA   Max.   : NA         Max.   : NA   Max.   : NA   Max.   : NA  
#>     evenness    jaccard_dis    nestedness     turnover     first_year 
#>  Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA   Min.   : NA  
#>  1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA   1st Qu.: NA  
#>  Median : NA   Median : NA   Median : NA   Median : NA   Median : NA  
#>  Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA   Max.   : NA  
#>    last_year        span    
#>  Min.   : NA   Min.   : NA  
#>  1st Qu.: NA   1st Qu.: NA  
#>  Median : NA   Median : NA  
#>  Mean   :NaN   Mean   :NaN  
#>  3rd Qu.: NA   3rd Qu.: NA  
#>  Max.   : NA   Max.   : NA

Tremendous trends in abundance:

high_abundance_slope <- slope_df %>%
  arrange(desc(abs(total_abundance)))
high_abundance_slope
#> # A tibble: 2,678 × 22
#>    siteid total_abundance log_total_abundance species_nb log_species_nb
#>    <chr>            <dbl>               <dbl>      <dbl>          <dbl>
#>  1 S11735           -756.              -0.971      0.245          1.31 
#>  2 S6466            -472.              -7.44       0.309          1.38 
#>  3 S5180             398.               3.91       0.207          0.822
#>  4 S5812            -368.               1.23       0.393          1.56 
#>  5 S3863             330.              12.6        0.206          1.53 
#>  6 S1787             277.              29.3        0.337          2.69 
#>  7 S8577             256.               4.78       0.168          0.800
#>  8 S4501             255.              29.4        0.582         10.8  
#>  9 S8158            -249.              -6.37       0.108          0.562
#> 10 S6074             245.               2.37       0.198          0.908
#> # … with 2,668 more rows, and 17 more variables: chao_richness <dbl>,
#> #   chao_shannon <dbl>, chao_simpson <dbl>, chao_evenness <dbl>, jaccard <dbl>,
#> #   horn <dbl>, chao <dbl>, hillebrand <dbl>, total <dbl>, appearance <dbl>,
#> #   disappearance <dbl>, evenness <dbl>, shannon <dbl>, simpson <dbl>,
#> #   jaccard_dis <dbl>, nestedness <dbl>, turnover <dbl>
site_desc_slope_total_abundance <-
  high_abundance_slope %>%
  .[["siteid"]]

The most tremendous temporal trends in abundance are located in
Illinois and France

ti <- filtered_dataset$location %>%
  filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 2154)

sp::plot(ne_countries())
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
  filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
  group_by(siteid) %>%
  summarise(across(c("total_abundance", "species_nb"), list(median = median, mean = mean)))
#> # A tibble: 20 × 5
#>    siteid total_abundance_me… total_abundance_… species_nb_medi… species_nb_mean
#>  * <chr>                <dbl>             <dbl>            <dbl>           <dbl>
#>  1 S11735               5948             17416.             18             18.7 
#>  2 S1507                  86               698.              2             10.4 
#>  3 S1787                2790.             2498.             16.5           15.7 
#>  4 S2595                 882              1306.             15             15.9 
#>  5 S3445                1791              2112.              5.5            5.31
#>  6 S3863                1827              3353.             13             13.6 
#>  7 S3883                 755              1842.             10              9.46
#>  8 S4325                 644.              899.              5.5            5.58
#>  9 S4501                2111              1924.             16             14.9 
#> 10 S4755                2517              2335.             20             20.1 
#> 11 S5180                4335              8431.             23             23.3 
#> 12 S5392                1416              4011.             19.5           18.9 
#> 13 S5812                4968             12216.             24.5           25.2 
#> 14 S5987                1470              1630.             18             14.2 
#> 15 S6074                5908.             9804.             23.5           23.1 
#> 16 S6466                3166              6465.             27             25.7 
#> 17 S8158                3610              5140.             24.5           24.7 
#> 18 S8577                4474              7435.             28             26.4 
#> 19 S8852                4126              5016.             25             22.8 
#> 20 S9031                3078             11660.             22             22.2

p_high_slope_abun <- map(site_desc_slope_total_abundance[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, total_abundance),
    y = "total_abundance", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_high_slope_abun, ncol = 3)

s11735_com <- analysis_dataset %>%
  filter(siteid == site_desc_slope_total_abundance[1]) %>%
  arrange(year) %>%
  select(siteid, year, total_abundance)
s11735_com
#> # A tibble: 25 × 3
#>    siteid  year total_abundance
#>    <chr>  <dbl>           <dbl>
#>  1 S11735  1993              66
#>  2 S11735  1994            1366
#>  3 S11735  1995          109088
#>  4 S11735  1996            9997
#>  5 S11735  1997             512
#>  6 S11735  1998            5948
#>  7 S11735  1999           85163
#>  8 S11735  2000           10785
#>  9 S11735  2001           21720
#> 10 S11735  2002           11719
#> # … with 15 more rows
plot_community_data(s11735_com, y = "total_abundance", x = "year",
  smoothing_method = "lm")

mean(s11735_com$total_abundance) -  sd(s11735_com$total_abundance)
#> [1] -10662.47
mean(s11735_com$total_abundance) +  5 * sd(s11735_com$total_abundance)
#> [1] 157805.5

Petit poissons Notropis atherinoides and Notropis volucellus have several tremendous abundance. Since those bursts happened four times, there have little probability to be error in the data.

s11735_pop <- filtered_dataset$measurement %>%
  filter(siteid == site_desc_slope_total_abundance[1]) %>%
  arrange(year) %>%
  select(siteid, year, species, abundance)

s11735_pop %>%
  group_by(species) %>%
  summarise(abun = sum(abundance)) %>%
  arrange(desc(abun))
#> # A tibble: 52 × 2
#>    species                   abun
#>    <chr>                    <dbl>
#>  1 Notropis atherinoides   418706
#>  2 Notropis volucellus       8400
#>  3 Cyprinella spiloptera     1949
#>  4 Lepomis macrochirus       1659
#>  5 Dorosoma cepedianum        651
#>  6 Aplodinotus grunniens      614
#>  7 Morone chrysops            540
#>  8 Pimephales vigilax         527
#>  9 Cyprinus carpio            427
#> 10 Macrhybopsis aestivalis    386
#> # … with 42 more rows

plot_temporal_population(s11735_pop) +
  theme(legend.position = "bottom")
#> Warning: Removed 52 rows containing missing values (position_stack).

Second highest temporal trends in abundance:

s6466_com <- analysis_dataset %>%
  filter(siteid == site_desc_slope_total_abundance[2]) %>%
  arrange(year) %>%
  select(siteid, year, total_abundance)
s6466_com
#> # A tibble: 11 × 3
#>    siteid  year total_abundance
#>    <chr>  <dbl>           <dbl>
#>  1 S6466   1981            7844
#>  2 S6466   1983           11938
#>  3 S6466   1987            9521
#>  4 S6466   1988           23997
#>  5 S6466   1989            1558
#>  6 S6466   1990            1674
#>  7 S6466   1991            3166
#>  8 S6466   1992             862
#>  9 S6466   1997            2687
#> 10 S6466   1999            1558
#> 11 S6466   2000            6314
plot_community_data(s6466_com, y = "total_abundance", x = "year",
  smoothing_method = "lm")

s6466_pop <- filtered_dataset$measurement %>%
  filter(siteid == site_desc_slope_total_abundance[2]) %>%
  arrange(year) %>%
  select(siteid, year, species, abundance)

s6466_pop %>%
  group_by(species) %>%
  summarise(abun = sum(abundance)) %>%
  arrange(desc(abun))
#> # A tibble: 55 × 2
#>    species                abun
#>    <chr>                 <dbl>
#>  1 Dorosoma cepedianum   43869
#>  2 Aplodinotus grunniens  8914
#>  3 Notropis wickliffi     8637
#>  4 Alosa chrysochloris    1812
#>  5 Morone chrysops        1538
#>  6 Cyprinus carpio        1272
#>  7 Ictalurus punctatus    1192
#>  8 Notropis atherinoides   745
#>  9 Ictiobus bubalus        559
#> 10 Sander canadensis       465
#> # … with 45 more rows

plot_temporal_population(s6466_pop) +
  theme(legend.position = "bottom")
#> Warning: Removed 495 rows containing missing values (position_stack).

3 Relationships among temporal trends

3.1 Globally

cor_slope <- slope_df %>%
  select(-siteid) %>%
  cor(., method = "spearman")
cor_slope
#>                     total_abundance log_total_abundance   species_nb
#> total_abundance         1.000000000         0.830317697  0.277156569
#> log_total_abundance     0.830317697         1.000000000  0.335207093
#> species_nb              0.277156569         0.335207093  1.000000000
#> log_species_nb          0.250125492         0.341531037  0.938150124
#> chao_richness           0.016699737         0.078327544  0.815025195
#> chao_shannon           -0.060971568        -0.004315524  0.598119627
#> chao_simpson           -0.068607673        -0.030703462  0.486160737
#> chao_evenness          -0.046333174         0.001842386  0.158532100
#> jaccard                 0.029465496         0.058720565 -0.023342671
#> horn                    0.032107656         0.065162610 -0.008994887
#> chao                    0.049896277         0.074583303  0.068878719
#> hillebrand              0.049560233         0.095796198 -0.025148726
#> total                   0.020626140         0.011421991  0.205181084
#> appearance              0.161650857         0.192922044  0.687912675
#> disappearance          -0.161078700        -0.222194738 -0.555815437
#> evenness               -0.135468004        -0.106903064  0.191546294
#> shannon                -0.030652568         0.031284774  0.595797755
#> simpson                -0.036762181         0.012861550  0.475730622
#> jaccard_dis            -0.028841716        -0.058365304  0.023323286
#> nestedness              0.004895141        -0.007242935  0.105858491
#> turnover               -0.022365349        -0.038795778 -0.052263750
#>                     log_species_nb chao_richness chao_shannon chao_simpson
#> total_abundance        0.250125492    0.01669974 -0.060971568  -0.06860767
#> log_total_abundance    0.341531037    0.07832754 -0.004315524  -0.03070346
#> species_nb             0.938150124    0.81502519  0.598119627   0.48616074
#> log_species_nb         1.000000000    0.79177426  0.652142713   0.55798584
#> chao_richness          0.791774260    1.00000000  0.715598168   0.58506236
#> chao_shannon           0.652142713    0.71559817  1.000000000   0.94429606
#> chao_simpson           0.557985838    0.58506236  0.944296056   1.00000000
#> chao_evenness          0.238494063    0.17861100  0.627245251   0.67332684
#> jaccard               -0.022860296   -0.03144932 -0.016411777  -0.02714028
#> horn                  -0.002333034   -0.01954284 -0.001651297  -0.01305221
#> chao                   0.078669559    0.05198781  0.064378306   0.04836406
#> hillebrand            -0.022070175   -0.03353840 -0.068107966  -0.08324337
#> total                  0.197786137    0.18733040  0.126450573   0.11238128
#> appearance             0.691348158    0.58925101  0.453612568   0.38451499
#> disappearance         -0.577843872   -0.45943156 -0.378710413  -0.30870107
#> evenness               0.289979041    0.32190903  0.722795822   0.77814032
#> shannon                0.652755093    0.67989789  0.944628419   0.91122873
#> simpson                0.555069101    0.55478316  0.873409696   0.89103134
#> jaccard_dis            0.022736545    0.03132624  0.016070172   0.02665117
#> nestedness             0.097239453    0.10838694  0.094533541   0.07288059
#> turnover              -0.048541424   -0.05721969 -0.064561731  -0.03107103
#>                     chao_evenness      jaccard         horn        chao
#> total_abundance      -0.046333174  0.029465496  0.032107656  0.04989628
#> log_total_abundance   0.001842386  0.058720565  0.065162610  0.07458330
#> species_nb            0.158532100 -0.023342671 -0.008994887  0.06887872
#> log_species_nb        0.238494063 -0.022860296 -0.002333034  0.07866956
#> chao_richness         0.178610995 -0.031449324 -0.019542841  0.05198781
#> chao_shannon          0.627245251 -0.016411777 -0.001651297  0.06437831
#> chao_simpson          0.673326844 -0.027140280 -0.013052206  0.04836406
#> chao_evenness         1.000000000  0.063506484  0.074803390  0.06164999
#> jaccard               0.063506484  1.000000000  0.989134887  0.69772222
#> horn                  0.074803390  0.989134887  1.000000000  0.73627697
#> chao                  0.061649988  0.697722217  0.736276969  1.00000000
#> hillebrand            0.014641849  0.332833278  0.347381250  0.24936211
#> total                -0.042854829 -0.933845455 -0.920529480 -0.59959876
#> appearance            0.090571956 -0.547707104 -0.527710202 -0.28291096
#> disappearance        -0.143240688 -0.642940230 -0.662605243 -0.52031520
#> evenness              0.712515794  0.047910610  0.058183834  0.06970014
#> shannon               0.565547657 -0.005489195  0.009607005  0.07637940
#> simpson               0.595960699 -0.002340240  0.012570159  0.06921570
#> jaccard_dis          -0.064111974 -0.999857414 -0.989090989 -0.69677143
#> nestedness           -0.008434623 -0.405843668 -0.355504999 -0.33184081
#> turnover             -0.037696468 -0.492788826 -0.522609730 -0.27616912
#>                      hillebrand       total  appearance disappearance
#> total_abundance      0.04956023  0.02062614  0.16165086   -0.16107870
#> log_total_abundance  0.09579620  0.01142199  0.19292204   -0.22219474
#> species_nb          -0.02514873  0.20518108  0.68791268   -0.55581544
#> log_species_nb      -0.02207017  0.19778614  0.69134816   -0.57784387
#> chao_richness       -0.03353840  0.18733040  0.58925101   -0.45943156
#> chao_shannon        -0.06810797  0.12645057  0.45361257   -0.37871041
#> chao_simpson        -0.08324337  0.11238128  0.38451499   -0.30870107
#> chao_evenness        0.01464185 -0.04285483  0.09057196   -0.14324069
#> jaccard              0.33283328 -0.93384545 -0.54770710   -0.64294023
#> horn                 0.34738125 -0.92052948 -0.52771020   -0.66260524
#> chao                 0.24936211 -0.59959876 -0.28291096   -0.52031520
#> hillebrand           1.00000000 -0.31748997 -0.21212267   -0.20672017
#> total               -0.31748997  1.00000000  0.72528811    0.48889404
#> appearance          -0.21212267  0.72528811  1.00000000   -0.11647853
#> disappearance       -0.20672017  0.48889404 -0.11647853    1.00000000
#> evenness            -0.05321078 -0.01433782  0.13964393   -0.17573304
#> shannon             -0.07563072  0.11786971  0.44713913   -0.38733767
#> simpson             -0.07581108  0.08965621  0.36745811   -0.32242970
#> jaccard_dis         -0.33261575  0.93400650  0.54788050    0.64308627
#> nestedness          -0.07413472  0.38337388  0.24465363    0.08227401
#> turnover            -0.21287472  0.47223951  0.30131921    0.50349434
#>                         evenness      shannon      simpson  jaccard_dis
#> total_abundance     -0.135468004 -0.030652568 -0.036762181 -0.028841716
#> log_total_abundance -0.106903064  0.031284774  0.012861550 -0.058365304
#> species_nb           0.191546294  0.595797755  0.475730622  0.023323286
#> log_species_nb       0.289979041  0.652755093  0.555069101  0.022736545
#> chao_richness        0.321909029  0.679897895  0.554783164  0.031326242
#> chao_shannon         0.722795822  0.944628419  0.873409696  0.016070172
#> chao_simpson         0.778140322  0.911228727  0.891031340  0.026651172
#> chao_evenness        0.712515794  0.565547657  0.595960699 -0.064111974
#> jaccard              0.047910610 -0.005489195 -0.002340240 -0.999857414
#> horn                 0.058183834  0.009607005  0.012570159 -0.989090989
#> chao                 0.069700135  0.076379395  0.069215704 -0.696771428
#> hillebrand          -0.053210777 -0.075630723 -0.075811084 -0.332615747
#> total               -0.014337815  0.117869711  0.089656212  0.934006505
#> appearance           0.139643932  0.447139132  0.367458109  0.547880499
#> disappearance       -0.175733040 -0.387337667 -0.322429697  0.643086266
#> evenness             1.000000000  0.784047749  0.866977324 -0.048397766
#> shannon              0.784047749  1.000000000  0.962316417  0.005099809
#> simpson              0.866977324  0.962316417  1.000000000  0.001943911
#> jaccard_dis         -0.048397766  0.005099809  0.001943911  1.000000000
#> nestedness           0.009438247  0.090073607  0.065840971  0.405983665
#> turnover            -0.053642052 -0.071596893 -0.052102504  0.492799673
#>                       nestedness    turnover
#> total_abundance      0.004895141 -0.02236535
#> log_total_abundance -0.007242935 -0.03879578
#> species_nb           0.105858491 -0.05226375
#> log_species_nb       0.097239453 -0.04854142
#> chao_richness        0.108386938 -0.05721969
#> chao_shannon         0.094533541 -0.06456173
#> chao_simpson         0.072880587 -0.03107103
#> chao_evenness       -0.008434623 -0.03769647
#> jaccard             -0.405843668 -0.49278883
#> horn                -0.355504999 -0.52260973
#> chao                -0.331840814 -0.27616912
#> hillebrand          -0.074134723 -0.21287472
#> total                0.383373875  0.47223951
#> appearance           0.244653628  0.30131921
#> disappearance        0.082274006  0.50349434
#> evenness             0.009438247 -0.05364205
#> shannon              0.090073607 -0.07159689
#> simpson              0.065840971 -0.05210250
#> jaccard_dis          0.405983665  0.49279967
#> nestedness           1.000000000 -0.47739313
#> turnover            -0.477393127  1.00000000

png(here("doc", "fig", "p_cor_slope.png"))
corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
      "total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)
dev.off()
#> png 
#>   2

png(here("doc", "fig", "p_cor_slope_tot.png"))
corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)
dev.off()
#> png 
#>   2

corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
      "total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)


corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)

ti <- expand.grid(
  resp1 = unique(slope$response),
  resp2 = unique(slope$response)
  ) %>%
  filter(resp2 != resp1) %>%
  filter(
    resp1 %in% c("chao_richness", "species_nb", "log_species_nb",
      "total_abundance", "log_total_abundance")) %>%
  mutate_all(as.character) %>%
  arrange(resp1)

test <- map2(ti$resp1, ti$resp2,
  function(x, y) {
    bi <- slope %>%
      filter(response %in% c(x, y)) %>%
      select(siteid, response, linear_slope) %>%
      pivot_wider(names_from = "response", values_from = "linear_slope")

    return(bi)
  }
)

p_trends_trends <- map(test, function(x) {
  l <- colnames(x)
  x %>%
    ggplot(aes(x = !!sym(l[2]), y = !!sym(l[3]))) +
    geom_point() +
    geom_smooth(method = "loess") +
    labs(
      x = get_var_replacement()[l[2]],
      y = get_var_replacement()[l[3]]
    )
})
names(p_trends_trends) <- map_chr(test, ~colnames(.x)[2])

3.2 Total abundance

plot_grid(
  plotlist = p_trends_trends[names(p_trends_trends) %in% "total_abundance"],
  ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

3.3 Species richness

plot_grid(
  plotlist = p_trends_trends[names(p_trends_trends) %in% "species_nb"],
  ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

3.4 Log species richness

plot_grid(
  plotlist = p_trends_trends[names(p_trends_trends) %in% "log_species_nb"],
  ncol = 3
)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

Trends of total abundance are linked to the one of species richness but not to shannon and simpson, suggesting that the trends of abundance are more link to trends of rare species. Logic: gain and loss of species should concern rare species.
Species richness decrease linked to disappearance, but not increase
Species richness increase linked to appearance, but not decrease

3.5 Turnover vs nestedness:

test <- slope_df %>%
  mutate(
    nes_sup_tu = nestedness > turnover,
    nes_inf_tu = nestedness < turnover,
    rich_inc = species_nb >= 0
  )
sum(test$nes_sup_tu)
#> [1] 1224
sum(test$nes_inf_tu)
#> [1] 1318
test %>%
  tabyl(rich_inc, nes_sup_tu) %>%
  adorn_totals("row") %>%
  adorn_title(placement = "top")
#>           nes_sup_tu     
#>  rich_inc      FALSE TRUE
#>     FALSE        733  499
#>      TRUE        721  725
#>     Total       1454 1224

median_site_rich <- analysis_dataset %>%
  select(c("siteid", starts_with("chao_"), "total_abundance")) %>%
  group_by(siteid) %>%
  summarise(across(where(is.double), median)) %>%
  rename_with(~paste0("med_", .x), where(is.double))

slope_df_med <- slope_df %>%
  left_join(median_site_rich, by = "siteid") %>%
  mutate(
    nes_sup_tu = nestedness > turnover,
    nes_inf_tu = nestedness < turnover,
    rich_inc = species_nb >= 0
  )

slope_df_med %>%
  ggplot(aes(
      x = log_species_nb,
      y = jaccard_dis,
      shape = nes_sup_tu,
      color = nes_sup_tu
      )) +
  geom_point() +
  geom_smooth(method = "gam")
#> `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

p_baselga_jaccard <- map(c("nestedness", "turnover"),
  ~slope_df_med %>%
    ggplot(aes(
        x = jaccard_dis,
        y = !!sym(.x),
      color = log(med_chao_richness))) +
    geom_point() +
    scale_color_viridis() +
    geom_abline(intercept = 0, slope = 1) +
    geom_smooth(method = "loess")
)
p_baselga_jaccard[[3]] <-
  slope_df_med %>%
  ggplot(aes(
      x = jaccard_dis,
      y = nestedness - turnover,
      color = log(med_chao_richness))) +
  geom_point() +
  scale_color_viridis() +
  geom_abline(intercept = 0, slope = 1) +
  geom_smooth(method = "loess")

plot_grid(plotlist = p_baselga_jaccard, ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

slope_df_med %>%
  ggplot(aes(x = med_chao_richness, y = jaccard_dis)) +
  geom_point()

4 Look at strong turnover trends

Check with Jaccard:

high_jaccard_slope <-
  slope_df %>%
  select(siteid, total_abundance, species_nb, jaccard) %>%
  arrange(desc(jaccard))
high_jaccard_slope
#> # A tibble: 2,678 × 4
#>    siteid total_abundance species_nb jaccard
#>    <chr>            <dbl>      <dbl>   <dbl>
#>  1 S450          120.        0.870    0.0673
#>  2 S6544           0.257    -0.0769   0.0559
#>  3 S9623          -3.77     -0.0266   0.0556
#>  4 S3217           0.0387   -0.0473   0.0496
#>  5 S6550           1.92      0.116    0.0490
#>  6 S4875          28.3       0.145    0.0461
#>  7 S2831           6.25     -0.0455   0.0420
#>  8 S9060          -0.817    -0.0303   0.0417
#>  9 S5514           0.690    -0.0576   0.0403
#> 10 S2742           1.08     -0.00909  0.0394
#> # … with 2,668 more rows
site_desc_jaccard_slope <-
  high_jaccard_slope %>%
  .[["siteid"]]

ti <- filtered_dataset$location %>%
  filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 2154)

sp::plot(ne_countries())
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
  filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
  group_by(siteid) %>%
  summarise(across(c("total_abundance", "species_nb"),
      list(median = median, mean = mean))) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
  scroll_box(width = "100%", height = "400px")

siteid	total_abundance_median	total_abundance_mean	species_nb_median	species_nb_mean
S11001	29.5700	33.575000	4.0	4.181818
S1574	323.0000	382.245454	4.0	3.272727
S2430	40.6520	39.779385	1.0	1.461539
S2742	41.9000	57.090909	2.0	2.000000
S2831	170.5000	180.154546	1.0	1.272727
S3044	40.5000	50.280000	1.0	1.100000
S3217	7.7000	8.793333	3.0	2.800000
S450	877.5000	744.666667	12.5	9.166667
S4875	290.5000	258.900000	13.0	11.200000
S5046	546.5000	596.100000	15.5	11.200000
S5514	25.9000	23.580000	1.5	1.700000
S560	44.1500	45.183333	3.0	3.416667
S6544	15.0000	18.061539	1.0	1.307692
S6550	6.6500	14.557143	2.0	2.214286
S6783	10.1000	12.246667	1.0	1.200000
S6936	22.8000	26.625000	1.0	1.333333
S7806	10.2000	10.950000	1.0	1.222222
S8447	10.7000	11.000000	4.0	3.750000
S9060	21.8185	21.676300	1.5	1.500000
S9623	61.0370	54.827200	2.5	2.500000

p_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
  ~plot_community_data(
    analysis_dataset %>%
    filter(siteid == .x) %>%
    select(siteid, year, jaccard),
    y = "jaccard", x = "year",
  smoothing_method = "lm"))
plot_grid(plotlist = p_high_jaccard_slope, ncol = 3)
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).

tar_load(vegdist_turnover_c)

Hum, do not understand why community are more similar at the end of the time-steps

p_pop_high_jaccard <-
  map(site_desc_jaccard_slope[1:20],
    ~plot_temporal_population(
      filtered_dataset$measurement %>%
        filter(siteid == .x)
      ) +
    theme(legend.position = "none")
    )
plot_grid(plotlist = p_pop_high_jaccard, ncol = 3)
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 77 rows containing missing values (position_stack).
#> Warning: Removed 8 rows containing missing values (position_stack).
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 3 rows containing missing values (position_stack).
#> Warning: Removed 8 rows containing missing values (position_stack).
#> Warning: Removed 60 rows containing missing values (position_stack).
#> Warning: Removed 2 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 26 rows containing missing values (position_stack).

Lets see communities that becomes more dissimilar:

p_jaccard_low_jaccard <-
  map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
    ~plot_community_data(
      analysis_dataset %>%
        filter(siteid == .x) %>%
        select(siteid, year, jaccard),
      y = "jaccard", x = "year",
      smoothing_method = "lm"))
plot_grid(plotlist = p_jaccard_low_jaccard, ncol = 3)
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).

p_pop_low_jaccard <-
  map(site_desc_jaccard_slope[
    (length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
    ~plot_temporal_population(
      filtered_dataset$measurement %>%
        filter(siteid == .x)
      ) +
    theme(legend.position = "none")
    )
plot_grid(plotlist = p_pop_low_jaccard, ncol = 3)
#> Warning: Removed 60 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 312 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 9 rows containing missing values (position_stack).
#> Warning: Removed 30 rows containing missing values (position_stack).
#> Warning: Removed 63 rows containing missing values (position_stack).
#> Warning: Removed 28 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 60 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 12 rows containing missing values (position_stack).
#> Warning: Removed 27 rows containing missing values (position_stack).
#> Warning: Removed 126 rows containing missing values (position_stack).

5 Comparison rigal, lm

rigal_trends_df <-  map2_dfr(
  rigal_trends, names(rigal_trends),
  ~.x %>% mutate(variable = .y)
  ) %>%
  select(variable, siteid, linear_slope) %>%
  pivot_wider(names_from = "variable", values_from = "linear_slope") %>%
  mutate(type = "rigal")

slope_comp <- list()
slope_comp$rigal <- rigal_trends_df[order(rigal_trends_df$siteid), ]
slope_comp$simple_lm <- slope_df[order(slope_df$siteid), ]
stopifnot(all(slope_comp$rigal$siteid == slope_comp$simple_lm$siteid))
slope_comp$diff <- tibble(
  siteid = slope_comp$rigal$siteid)
for (i in var_temporal_trends) {
  slope_comp$diff[[i]] <- slope_comp$rigal[[i]] - slope_comp$simple_lm[[i]]
}

All good !

old_par <- par()
par(mfrow = c(3, 2))
for (i in var_temporal_trends) {
  plot(
    slope_comp$rigal[[i]] ~
      slope_comp$simple_lm[[i]],
    ylab = "Rigal (polynomial degree 2)",
    xlab = "Simple LM"
  )
  abline(0, 1)
  title(i)
}

par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé

# Compare Rigal and simple lm for few sites
plot_rigal_raw_data <- function(
  site = NULL,
  response = NULL,
  dataset = NULL,
  rigal_trends = NULL
  ) {

  raw_data <- dataset %>%
    filter(siteid == site)
  rigal_resp <- rigal_trends[[response]]
  rigal_resp_site <- rigal_resp[rigal_resp$siteid == site, ]
  plot(raw_data[[response]]~raw_data[["year"]], ylab = response, xlab = "Year")
  abline(rigal_resp_site$intercept, rigal_resp_site$linear_slope, col
  = "red")
  abline(lm(raw_data[[response]]~raw_data[["year"]]), col = "green")
  title(main = site)
}

plot_rigal_raw_data(site = "S2672", response = "log_total_abundance", dataset =
analysis_dataset, rigal_trends = rigal_trends)

slope_comp$diff <- slope_comp$diff %>%
  arrange(desc(abs(log_total_abundance)))

old_par <- par()
par(mfrow = c(3, 2))
for (i in slope_comp$diff[1:6,]$siteid) {
  plot_rigal_raw_data(
    site = i,
    response = "log_total_abundance",
    dataset = analysis_dataset,
    rigal_trends = rigal_trends
  )
}

par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé

6 Classification of temporal trends

Classification:

classification_station <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = shape_class)) +
  geom_bar() +
  labs(x = "Trajectory classification", y = "Nb station",
    title = get_var_replacement()[.y]) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
)


plot_grid(plotlist = classification_station)

p_classif_abun_rich <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "abundance|richness|species")]])
save_plot(
  here("doc", "fig", "classif_trends_abun_rich.png"),
  p_classif_abun_rich
)
p_classif_hill_evenness <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[str_detect(names(classification_station),
  "shannon|simpson|evennes")]])
save_plot(
  here("doc", "fig", "classif_trends_hill_evenness.png"),
  p_classif_hill_evenness
)
p_classif_turnover <- plot_grid(
  plotlist = classification_station[
    names(classification_station)[names(classification_station) %in%
      c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
        "disappearance")]])
save_plot(
  here("doc", "fig", "classif_trends_classif_hill_evenness.png"),
  p_classif_turnover
)

ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "direction"),
  .id = "response"
)
ti %>%
  filter(direction != "Total") %>%
  select(response, direction, percent) %>%
  mutate(response = str_replace_all(response, get_var_replacement())) %>%
  pivot_wider(names_from = "direction", values_from = "percent") %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "bordered", "hover"))

response	stable	increase	decrease
Total abundance	75.8%	14.3%	9.9%
Log Total Turnover (jaccard) (codyn) abundance	73.3%	16.5%	10.2%
Species richness	77.7%	13.2%	7.4%
Log species richness	76.2%	13.2%	6.6%
Chao species richness	79.8%	10.6%	8.2%
Chao Shannon	76.5%	12.1%	7.7%
Chao Simpson	78.4%	10.5%	7.4%
Chao Evenness	82.8%	7.5%	5.9%
Jaccard (binary, similarity)	76.8%	2.8%	17.5%
Horn (binary, similarity)	77.3%	2.7%	17.2%
Chao (binary, similarity)	71.5%	2.3%	8.3%
SER_a (rel abundance)	75.3%	4.1%	18.6%
Total Turnover (jaccard) (codyn)	73.7%	18.3%	3.0%
Appearance	70.2%	15.6%	3.8%
disAppearance	64.0%	9.5%	3.0%
Evenness	78.8%	9.4%	8.1%
Shannon	75.7%	12.2%	8.3%
Simpson	77.9%	10.8%	7.7%
Jaccard (binary, dissimilarity)	74.7%	17.4%	2.8%
Nestedness (jaccard)	80.3%	9.3%	5.3%
Turnover (jaccard)	57.4%	8.3%	1.7%

Linear trends:

#3.3   Log-linear model:logYi=α+βXi+εiIn

#the log-linear model, the literal interpretation of the estimated
#coefficientˆβis that a one-unitincrease inXwill produce an expected increase in
#logYofˆβunits. 
#In terms ofYitself, this
#meansthat  the  expected  value  ofYis  multiplied  byeˆβ.   So  in  terms  of
#effects  of  changes  inXonY(unlogged):•Each 1-unit increase inXmultiplies the
#expected value ofYbyeˆβ.•To compute the effects onYof another change inXthan an
#increase of one unit, call thischangec, we need to includecin the exponent.  The
#effect of ac-unit increase inXis tomultiply the expected value ofYbyecˆβ. So the
#effect for a 5-unit increase inXwould bee5ˆβ.•For small values ofˆβ,
#approximatelyeˆβ≈1+ˆβ. We can use this for the following approxima-tion for a
#quick interpretation of the coefficients:  100·ˆβis the expected percentage
#changeinYfor a unit increase inX.  For instance forˆβ=.06,e.06≈1.06, so a 1-unit
#change inXcorresponds to (approximately) an expected increase inYof 6%

#https://data.library.virginia.edu/interpreting-log-transformations-in-a-linear-model/
hist_linear_plot <- map2(rigal_trends, names(rigal_trends),
  ~.x %>%
  ggplot(aes(x = linear_slope)) +
  geom_histogram() +
  labs(x = "Linear slope", y = "Number of station",
    title = get_var_replacement()[.y]) +
  ylim(c(0, 1000))
)
plot_grid(plotlist = hist_linear_plot)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).

type_response <- list(
  abun_rich = str_detect(names(classification_station),
    "abundance|richness|species"),
  hill_evenness = str_detect(names(classification_station),
    "shannon|simpson|evennes"),
  turnover = names(classification_station) %in%
    c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
      "disappearance")
)

map2(type_response, names(type_response),
  function(x, y) {
  p <- plot_grid(
    plotlist = hist_linear_plot[
      names(hist_linear_plot)[x]
      ]
  )
  save_plot(
    here("doc", "fig", paste0("hist_linear_slope_", y, ".png")),
    p
  )

  }
)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> $abun_rich
#> NULL
#> 
#> $hill_evenness
#> NULL
#> 
#> $turnover
#> NULL

7 Spatialization of temporal trends

rigal_trends_df_loc <-
  rigal_trends_df %>%
  left_join(
    select(filtered_dataset$location, siteid, ecoregion, country),
    by = c("siteid"))

No obvious geographic patterns

rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = log_species_nb, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = log_total_abundance, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = total, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = jaccard, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
  ggplot(aes(x = country, y = jaccard_dis, fill = ecoregion)) +
  geom_boxplot()


rigal_trends_df_loc %>%
    ggplot(aes(x = country, y = hillebrand, fill = ecoregion)) +
    geom_boxplot()

rigal_trends_df_loc <- filtered_dataset$location %>%
  left_join(rigal_trends_df, by = "siteid") %>%
  st_as_sf(coords = c("longitude", "latitude"),
  crs = 4326)

library(colorspace)
div_pal <- colorspace::sequential_hcl(11, "plasma")
world <- ne_countries(returnclass = "sf") 

bb <- st_bbox(rigal_trends_df_loc)

map_world_trends(df = rigal_trends_df_loc, y = "log_species_nb")

map_world_trends(df = rigal_trends_df_loc, y = "jaccard")

map_world_trends(df = rigal_trends_df_loc, y = "log_total_abundance")

rigal_trends_df %>%
  summarise(across(where(is.numeric), median, na.rm = T)) %>%
  unlist()
#>     total_abundance log_total_abundance          species_nb      log_species_nb 
#>        1.129827e-01        4.993106e-01        2.026907e-03        4.980022e-02 
#>       chao_richness        chao_shannon        chao_simpson       chao_evenness 
#>        2.926352e-04        7.389185e-04        7.189205e-04        0.000000e+00 
#>             jaccard                horn                chao          hillebrand 
#>       -3.750972e-03       -2.731569e-03       -1.598283e-17       -3.218861e-03 
#>               total          appearance       disappearance            evenness 
#>        3.225717e-03        9.691423e-04        0.000000e+00        0.000000e+00 
#>             shannon             simpson         jaccard_dis          nestedness 
#>        4.650060e-04        1.071122e-04        3.750972e-03        2.880739e-04 
#>            turnover 
#>        0.000000e+00

median(exp(rigal_trends_df$log_total_abundance) - 1)  * 100
#> [1] 64.75967
sd(exp(rigal_trends_df$log_total_abundance) - 1) * 100
#> [1] 3.425009e+24
median(exp(na.omit(rigal_trends_df$log_species_nb))  - 1) * 100
#> [1] 5.106305
sd(exp(na.omit(rigal_trends_df$log_species_nb)) - 1) * 100
#> [1] 1.019932e+12

8 Sensibility analysis of temporal trends

tar_load(c(rigal_trends_avg3y))
names(rigal_trends_avg3y) <- var_temporal_trends

tar_load(c(analysis_dataset_avg3y, measurement_avg3y, measurement))
tar_load(c(turnover_avg3y_c))
tar_load(c(vegdist_turnover_c))
tar_load(c(hillnb_avg3y))

slope_avg3y <- map_dfr(rigal_trends_avg3y[
  ! names(rigal_trends_avg3y) %in% c("jaccard_dis", "nestedness", "turnover")], ~.x %>%
  select(siteid, linear_slope),
.id = "response"
)

Number of sites that are in the different datasets:

u_site <- map(list(classic = slope, avg_3y = slope_avg3y), ~unique(.x$siteid))
sum(u_site[[1]] %in% u_site[[2]])
#> [1] 2678
sum(u_site[[2]] %in% u_site[[1]])
#> [1] 2678

comparison_baseline <- rbind(
  slope %>%
    mutate(type = "classic") %>%
    filter(siteid %in% u_site[["avg_3y"]]),
  slope_avg3y %>%
    mutate(type = "avg_3y")
  ) %>%
pivot_wider(names_from = type, values_from = linear_slope)

It is not too bad, but we can see some differences:
- slope negatively biaised for log_species_nb in 3 years baseline
- Bunch of trends in appearance/disappearance trends disappears with
  3 years baseline
- But note that in three years baseline, there are three points less data

comparison_baseline %>%
  ggplot(aes(x = classic, y = avg_3y)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  facet_wrap(~response, scale = "free")
#> Warning: Removed 8034 rows containing missing values (geom_point).

8.1 Reproducibility

Reproducibility receipt

## datetime
Sys.time()
#> [1] "2022-01-20 10:22:18 CST"

## repository
if(requireNamespace('git2r', quietly = TRUE)) {
  git2r::repository()
} else {
  c(
    system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
    system2("git", args = c("remote", "-v"), stdout = TRUE)
  )
}
#> Local:    main /home/alain/Documents/post-these/isu/RivFishTimeBiodiversityFacets
#> Head:     [3695d10] 2022-01-20: fix hillebrand: change entry matrix

## session info
sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 10 (buster)
#> 
#> Matrix products: default
#> BLAS:   /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRblas.so
#> LAPACK: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
#>  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
#>  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] parallel  stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] colorspace_2.0-0        betapart_1.5.4          untb_1.7-4             
#>  [4] iNEXT_2.0.20            factoextra_1.0.7        ade4_1.7-16            
#>  [7] ecocom_0.0.0.9000       inlabru_2.3.1           INLA_21.11.22          
#> [10] sp_1.4-5                foreach_1.5.1           Matrix_1.3-2           
#> [13] future_1.21.0           vegan_2.5-7             lattice_0.20-41        
#> [16] permute_0.9-5           codyn_2.0.5             janitor_2.1.0          
#> [19] viridis_0.5.1           viridisLite_0.3.0       cowplot_1.1.1          
#> [22] rnaturalearthdata_0.1.0 rnaturalearth_0.1.0     sf_1.0-4               
#> [25] rmarkdown_2.11          scales_1.1.1            kableExtra_1.3.1       
#> [28] here_1.0.1              lubridate_1.7.9.2       magrittr_2.0.1         
#> [31] forcats_0.5.1           stringr_1.4.0           dplyr_1.0.4            
#> [34] purrr_0.3.4             readr_2.1.1             tidyr_1.1.2            
#> [37] tibble_3.1.6            ggplot2_3.3.3           tidyverse_1.3.0        
#> [40] tarchetypes_0.3.2       targets_0.8.1           conflicted_1.1.0       
#> 
#> loaded via a namespace (and not attached):
#>   [1] snow_0.4-4         readxl_1.3.1       backports_1.2.1   
#>   [4] fastmatch_1.1-3    corrplot_0.84      plyr_1.8.6        
#>   [7] igraph_1.2.6       splines_4.0.5      gmp_0.6-2.1       
#>  [10] listenv_0.8.0      digest_0.6.27      htmltools_0.5.1.1 
#>  [13] fansi_0.5.0        memoise_2.0.0      cluster_2.1.1     
#>  [16] tzdb_0.2.0         openxlsx_4.2.3     globals_0.14.0    
#>  [19] modelr_0.1.8       prettyunits_1.1.1  rvest_0.3.6       
#>  [22] ggrepel_0.9.1      haven_2.3.1        xfun_0.28         
#>  [25] rgdal_1.5-28       callr_3.7.0        crayon_1.4.2      
#>  [28] jsonlite_1.7.2     ape_5.5            iterators_1.0.13  
#>  [31] glue_1.5.1         gtable_0.3.0       webshot_0.5.2     
#>  [34] car_3.0-10         abind_1.4-5        emo_0.0.0.9000    
#>  [37] DBI_1.1.1          rstatix_0.7.0      Rcpp_1.0.6        
#>  [40] progress_1.2.2     magic_1.5-9        units_0.6-7       
#>  [43] foreign_0.8-81     httr_1.4.2         RColorBrewer_1.1-2
#>  [46] wk_0.5.0           ellipsis_0.3.2     partitions_1.10-4 
#>  [49] pkgconfig_2.0.3    farver_2.0.3       sass_0.3.1        
#>  [52] dbplyr_2.1.0       utf8_1.2.2         polynom_1.4-0     
#>  [55] reshape2_1.4.4     tidyselect_1.1.1   labeling_0.4.2    
#>  [58] rlang_0.4.12       munsell_0.5.0      cellranger_1.1.0  
#>  [61] tools_4.0.5        cachem_1.0.4       cli_3.1.0         
#>  [64] generics_0.1.0     mathjaxr_1.4-0     broom_0.7.4       
#>  [67] geometry_0.4.5     evaluate_0.14      fastmap_1.1.0     
#>  [70] yaml_2.2.1         processx_3.5.2     knitr_1.36        
#>  [73] fs_1.5.1           zip_2.1.1          s2_1.0.7          
#>  [76] nlme_3.1-152       mime_0.11          xml2_1.3.2        
#>  [79] compiler_4.0.5     rstudioapi_0.13    curl_4.3.2        
#>  [82] e1071_1.7-4        ggsignif_0.6.1     reprex_1.0.0      
#>  [85] bslib_0.2.4        stringi_1.7.6      highr_0.9         
#>  [88] ps_1.6.0           Brobdingnag_1.2-6  rgeos_0.5-5       
#>  [91] classInt_0.4-3     vctrs_0.3.8        pillar_1.6.4      
#>  [94] lifecycle_1.0.1    furrr_0.2.2        jquerylib_0.1.3   
#>  [97] data.table_1.13.6  R6_2.5.1           bookdown_0.24     
#> [100] KernSmooth_2.23-18 gridExtra_2.3      rio_0.5.16        
#> [103] parallelly_1.23.0  codetools_0.2-18   rcdd_1.5          
#> [106] MASS_7.3-53.1      assertthat_0.2.1   picante_1.8.2     
#> [109] rprojroot_2.0.2    withr_2.4.3        doSNOW_1.0.19     
#> [112] mgcv_1.8-34        hms_1.1.1          grid_4.0.5        
#> [115] class_7.3-18       snakecase_0.11.0   carData_3.0-4     
#> [118] git2r_0.29.0       ggpubr_0.4.0       itertools_0.1-3

Temporal trends

Alain Danet

20 janvier, 2022

1 Summary

2 Check trends

3 Relationships among temporal trends

3.1 Globally

3.2 Total abundance

3.3 Species richness

3.4 Log species richness

3.5 Turnover vs nestedness:

4 Look at strong turnover trends

5 Comparison rigal, lm

6 Classification of temporal trends

7 Spatialization of temporal trends

8 Sensibility analysis of temporal trends

8.1 Reproducibility